#include<stdio.h>
#include<stdlib.h>
#include<limits.h>
#include<complex.h>
#include<fftw3.h>
#include<string.h>
#include<assert.h>
#include<math.h>
#include<time.h>

#define NENTRIES 64
#define MASK 0x3f
#define NENTRIES_CMP (NENTRIES+2)
double a_orig[NENTRIES];
double b_orig[NENTRIES];
double a[NENTRIES];
double b[NENTRIES];
double c1[NENTRIES];
double c2[NENTRIES];

void
cconv(double *a, double *b, double *c)
{
   double *d, *x, *y;
   double acc;
   int i,j;
   for (x=a, y=b, d=c,
        i=0; i<NENTRIES; i++) {
      acc = 0.;
      for (j=0; j<NENTRIES; j++) {
         acc += x[j] * y[(i-j)&MASK];
      }
      *d++ += acc;
   }
}

static int did_init=0;
fftw_plan pln_r2c_a;
fftw_plan pln_r2c_b;
fftw_plan pln_c2r_c;
double *inout_ac, *in_b;
fftw_complex *cmp_ac, *cmp_b;

void
cconv_fft(double *a, double *b, double *c)
{
   if (!did_init) {
      inout_ac = (double*) fftw_malloc(sizeof(double)*NENTRIES);
      in_b = (double*) fftw_malloc(sizeof(double)*NENTRIES);
      cmp_ac = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*NENTRIES_CMP);
      cmp_b = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*NENTRIES_CMP);
      assert(inout_ac && in_b && cmp_ac && cmp_b);
      pln_r2c_a = fftw_plan_dft_r2c_1d(NENTRIES, inout_ac, cmp_ac, FFTW_ESTIMATE);
      pln_r2c_b = fftw_plan_dft_r2c_1d(NENTRIES, in_b, cmp_b, FFTW_ESTIMATE);
      pln_c2r_c = fftw_plan_dft_c2r_1d(NENTRIES, cmp_ac, inout_ac, FFTW_ESTIMATE);
      did_init = 1;
   }
   /* copy in vectors */
   memcpy(inout_ac, a, sizeof(double)*NENTRIES);
   memcpy(in_b, b, sizeof(double)*NENTRIES);
   /* fft, mult in freq domain, inv-fft */
   fftw_execute(pln_r2c_a);
   fftw_execute(pln_r2c_b);
   int i;
   for (i=0; i<NENTRIES; i++) {
      cmp_ac[i] = cmp_ac[i] * cmp_b[i];
   }
   fftw_execute(pln_c2r_c);
   /* scale back down */
   for (i=0; i<NENTRIES; i++) {
      inout_ac[i] = inout_ac[i] / NENTRIES;
   }
   /* copy back result (in inout_ac) */
   memcpy(c, inout_ac, sizeof(double)*NENTRIES);
}

void
norm(double *v)
{
   double *cur;
   double len;
   int j;

   len = 0.;
   for (cur=v,
        j=0; j<NENTRIES; j++, cur++)
      len += (*cur) * (*cur);
   len = sqrt(len);
   for (cur=v,
        j=0; j<NENTRIES; j++, cur++)
      if (*cur != 0.)
         *cur /= len;
}

double
frand(void)
{
   return ((double)random())/(double)INT_MAX;
}

int main()
{
   srandom((unsigned)time(0));
   /* init a and b */
   int i;
   for (i=0; i<NENTRIES; i++) {
      a_orig[i] = frand(); //(double)i;
      b_orig[i] = frand(); //(double)i/2.;
      c1[i] = 0.;
      c2[i] = 0.;
   }
   for (i=0; i<NENTRIES; i++) {
      a[i] = a_orig[i];
      b[i] = b_orig[i];
   }
   //norm(a);
   //norm(b);
   cconv(a,b,c1);
   //norm(c1);
   /* repeat */
   for (i=0; i<NENTRIES; i++) {
      a[i] = a_orig[i];
      b[i] = b_orig[i];
   }
   //norm(a);
   //norm(b);
   cconv_fft(a,b,c2);
   //norm(c2);
   /* show results */
   for (i=0; i<NENTRIES; i++) {
      if (c1[i]!=c2[i]) 
         printf("%i %f %f (%f)\n", i, c1[i], c2[i], c2[i]-c1[i]);
   }
   return 0;
}
